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ABSTRACT 


We consider here the use of extra collecting horns in the focal plane of an antenna as a means 
of determining the Direction Of Arrival (DOA) of the signal impinging on it, provided it is within the 
antenna beam. Our analysis yields a relatively simple algorithm to extract the DOA from the horns’ 
outputs. We also develop here an algorithm which, in effect, measures the thermal noise of the horns’ 
signals and determines its effect on the uncertainty of the extracted DOA parameters. 

Both algorithms have been implemented in software and tested on simulated data. Based on 
these tests, we conclude that this is a viable approach to the DOA determination. 

Though the results obtained here are of general applicability, the particular motivation for the 
present work is their application to the pointing of a mechanically deformed antenna. It is anticipated 
that the pointing algorithm for a deformed antenna could be obtained as a small perturbation of the 
algorithm developed here for an undeformed antenna. In this context, it should be pointed out that, 
with a deformed antenna, the array of horns and its associated circuitry constitute the main part of the 
deformation-compensation system. In this case, the pointing system proposed here may be viewed as 
an additional task carried out by the deformation-compensation hardware. 
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I. 


INTRODUCTION 


Consider a circular-paraboloid reflector antenna in its reception mode. We have here a large 
"dish" which collects the energy arriving from a distant source and concentrates it (directly or via a 
Cassegrain system) onto the single collecting horn, which provides the antenna output. Conceptually, 
one can distinguish two basic planes in this system: the aperture plane, where the field has the 
uniform distribution of a plane wave, and the focal plane, where the field distribution is quite 
different, being characterized by a very pronounced peak at the location of the collecting horn. 

With a properly designed antenna, most of the energy incident on the antenna aperture is 
concentrated into this peak and thus fed to the output via the collecting horn. When such an antenna 
is subjected to a mechanical deformation, the field distribution in the focal plane will be diffuse so 
that the collecting horn will collect only part of the total energy incident on the antenna and we 
observe a signal loss. 

While the source of the mechanical deformation is of no concern to us here, we should point 
out a practically important scenario where such deformation has to be dealt with. We are referring to 
the situation where we try to extend upward the frequency coverage of an already-built functional 
antenna. In this case, the antenna whose mechanical deviation from the ideal is acceptable within its 
design frequency range will now be labeled "mechanically deformed" at the shorter wavelengths. 

One possible way of compensating for at least part of the deformation signal loss is to 
augment the single collecting horn with an array of additional horns closely packed around it and 
produce the overall system output as a weighted sum of all the horn outputs [1]. The main challenge 
in such a design is to determine the optimal combining weights and implement them in real time. 

When such a system is properly pointed to extract the maximum power from a desired source, 
it may be considered equivalent to a standard undeformed antenna of a somewhat smaller aperture, 
and one tends to presume that the deformation and the array of horns compensating for it are totally 
transparent to the user. This, however, is not the case since the behavior of the system when it is 
turned slightly away from this optimal direction is quite different from that of a standard antenna. For 
example, the mechanical deformation might be such that maximal output is produced when the 
antenna is pointed slightly off the actual direction of the source. This, then, raises a basic question: 
Given the direction of the desired source, how will such an antenna system be pointed to extract the 
maximal power from it? 

Obviously, the availability of the outputs of the individual horns of the collecting array 
provides the basis for solving this problem. We propose here a two-phase approach: In phase 1 we 
assume that the deformation is negligible and develop an algorithm which processes the outputs of the 
individual horns to yield the Direction Of Arrival (DOA) of a plane wave impinging on the antenna. 
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In phase 2 we consider the case of the actual deformed antenna as a small perturbation of the 
undeformed case and look for the modifications and additional parameters that have to be provided so 
that the DOA can still be extracted. In this report, we cover phase 1. 

Throughout this report we use the acronym DOA to represent the entity we want to extract, 
namely, the Direction Of Arrival of the signal. Strictly speaking, this is a precise characterization. 
However, there is a large amount of literature in which "extracting the DOA" has the additional 
implication of starting with a blank page, that is, we do not have even a coarse a priori estimate of 
the DOA. We want to stress at the outset that this is not the usage intended here. The algorithm 
developed in this report is applicable only when the antenna is already pointed at the vicinity of the 
source. Thus, this algorithm handles the task of fine pointing after coarse pointing has already been 
effected by other means. 

The general outline of the report is as follows: In Section II we develop the basic strategy 
underlying the DOA extraction. We show there that the DOA is determined by the gradient of the 
phase angle of the field entities across the antenna aperture. Since the aperture-plane entities are not 
directly available to us, we have to infer them from the focal-plane "measurements," namely, the 
outputs of the focal-plane horns. This is tackled in Section III. 

From this point on, our analysis is strongly influenced by the fact that the signals available to 
us are noisy. To (partially) overcome the effect of this noise, we have to increase the number of 
measurements and are thus naturally led to a least-squares formulation of the problem (number of 
measurements exceeds number of unknowns). It turns out that, in the least-squares context, it is more 
convenient to regard the DOA extraction as the determination of the orientation of a wavefront of the 
incident wave. The linear operator underlying this formulation is derived in Section IV. In Section V 
we study the structure of this linear operator using the powerful concepts of Singular Value 
Decomposition (SVD). We find that the representation of this operator yielded by our formulation is 
so simple that it allows us to obtain a closed-form solution of the least-squares problem (no matrix 
inversion; no SVD subroutines). 

At this point of the analysis we are in possession of an algorithm to obtain the best 
(least-squares) estimate of the DOA. The next three sections are devoted to the much more 
challenging task of determining the reliability of this estimate. We express this reliability as follows: 

A 

The estimated DOA is representable as a point (D) on the unit sphere. We look for the shape and size 
of the region around this point in which the point D representing the true DOA is guaranteed to fall 
with a prescribed probability. 

The first step towards this goal is the determination of the covariance matrix of the 
least-squares solution vector. In Section VI we express this covariance matrix in terms of the 
covariance matrix of the parameters at the input of the least-squares problem. Here we encounter a 
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major challenge: The input parameters of the least-squares problem are the phases of the 
aperture-plane field entities. But the noise is injected in the focal plane and the transformation 
between these two planes transforms the diagonal noise covariance matrix of the focal-plane into a 
nondiagonal noise covariance matrix in the aperture plane where the least-squares processing is 
applied. 

Section VII is devoted to the determination of this nondiagonal noise covariance matrix. The 
mathematical manipulations in this section are quite complex, but we have provided here a sufficiently 
detailed step-by-step derivation to compensate for that. 

In Section VIII we handle the extraction of the confidence region from the (now available) 
output covariance matrix. The final result here is an algorithm to determine the confidence region for 
any probability. 

The two algorithms developed in this report have been implemented in software and 
incorporated in a very flexible program which allows their testing with simulated or real data, with or 
without added noise, and with the easy selection of a large number of operational parameters. 
Preliminary results with simulated data indicate that these algorithms constitute a technically sound 
approach to the DOA determination. Some of these results are discussed in Section IX. 

In Section X we present a brief overview and a discussion of some worthwhile extensions of 
this work. 

Throughout this report, we use the formalism of dyadics to handle linear operators. When 
combined with the concepts of SVD, this provides a very powerful tool which has proved to be quite 
useful in mathematical manipulations as well as on the conceptual level. For those not familiar with 
this formalism, we provide a short but sufficient summary of the subject in Appendix A. 

Some of the detailed arguments in Section III require knowledge of the attenuation of the 
focal-plane field as a function of distance from the optical axis when the incident plane wave is off 
boresight. In Appendix B we derive this result. An incidental by-product of this derivation is the 
possibility that it may provide the basis for a DOA algorithm which is quite different from the one 
developed in this report. This is briefly discussed in this appendix. 

The number of variables in this work is quite large. To reduce the number of required 
symbols, we have made extensive use of the caret (*) to allow use of the same letter for two different 
but related entities. Even so, we find that some symbols have to be assigned different meanings in the 
different sections. The quite different contexts, however, should rule out any ambiguity. 
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n. 


THE BASIC STRATEGY 


The due to the DOA extraction is the very special variation of the field entities of a plane 
wave as functions of time and space. Specifically, for a plane wave of radial frequency o> (and 
corresponding wavelength A, ) , all components of the electric and magnetic fields vary as 

q i (k-r- o) t) (2.1) 

where r is the position vector of the point where the field is specified (throughout this report, we use 
overbars to denote vectors), t is the relevant time, and the vector k is the propagation vector, namely, 
a vector of magnitude k = 2%/X pointing in the direction of propagation of the plane wave. (Note: 
we adopt here the convention that the symbol of a vector entity without the overbar stands for the 
magnitude of that vector). 

Note that the dimensionless entity 


appearing in (2.1) is a phase angle. Thus, if we have sensors which continuously monitor a certain 
component of the field in various locations in space, then all of them will display a sinusoidal 
variation with time with the same radial frequency but with different phases, the phase of each sensor 
being given by the a corresponding to its position vector r (see (2.2)). The DOA we are after will be 
specified in terms of a unit vector ( n ) pointing at the source of the wave. From the definition 
of k , it is obvious that 

k - ~ (2 ir/A) n ^ 3) 

and our problem is essentially that of determining k . With that in mind, we return to (2.2) and 
consider the gradient of the phase angle a , namely, 

Va = k • Vr ( 2 . 4 ) 


But 


Vf = J 


( 2 . 5 ) 
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where J is the identity dyadic. (Throughout this report, we use double overbars to denote dyadics. 
See Appendix A for a short but sufficient summary of the subject of dyadics as well as a simple 
derivation of (2.5).) Hence, 


k = Va 


( 2 . 6 ) 


and, using (2.3), 



(2.7) 


We conclude that measuring the phase differences across an array of sensors monitoring the 
plane wave will yield the desired DOA. This, then, specifies the main points of our basic strategy: 

(1) We consider the focal-plane horns as sensors monitoring the focal-plane field. 

(2) We use the Fourier-transform relationship between the aperture-plane field and the 
focal-plane field to obtain from the horns’ outputs a field distribution across the 
aperture plane. 

(3) Applying (2.7) to the phases of the aperture-plane distribution, we obtain the desired 
DOA. 

In principle, (2.7) requires a three-dimensional phase distribution rather than the planar 
distribution obtained in the scheme outlined above. However, in our case, a planar distribution is 
sufficient, as we proceed to show now. 

Let us introduce an antenna-fixed Cartesian frame {Aj} j Ml as shown in Fig. 1. In this frame, 
all aperture-plane points are characterized by zero ii 3 components. This motivates us to split Va as 
follows: 


Va = (h 2 h 3 ) • (Va) + (J - h 3 h 3 ) • (Va) 


( 2 . 8 ) 


Let us denote 


I - h 2 h 2 ■ P 3 


(2.9) 
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P 3 is a projection operator which extracts from any vector its part orthogonal to h 2 . Thus, 


Va = (h^h^ • (Va) + P 3 • (Va) 

and the entity available to us from the aperture-plane distribution is notVa but 


( 2 . 10 ) 


P 3 • (Va) s V s a 


( 2 . 11 ) 


Such an entity is usually referred to as a "surface gradient." 

To see the effect of this shortcoming, let us apply P 3 to both sides of (2.7): 


P 


3 



Now we express n in terms of the h i frame 


n 


3 


= 52 n i h 

i - 1 


i 


( 2 . 12 ) 


(2.13) 


Hence, 


n = n 1 h 1 + n 2 h 2 = n. 


(2.14) 


and the final result is 


n s = Vo a 

s 2rt s 


(2.15) 


Note that 


n = n s + n 3 h^ 


(2.16) 
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with 


n 


3 


± 



(2.17) 


Thus, once we have computed n s , the location of the source is constrained to just two 
mathematically possible alternatives, one in front of the antenna ( n 3 > 0) and one behind 
it ( n 3 < 0) . Since the choice between these two is obvious, we conclude that the surface gradient 
is sufficient to determine the DOA in our case. 

In the next section, we take up step (2) of our basic strategy, namely, determining the 
unavailable aperture-plane distribution from the available focal-plane distribution. 
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m. EXTRACTING APERTURE-PLANE FIELD PARAMETERS FROM FOCAL-PLANE 


MEASUREMENTS 

Consider the field distributions in the two basic planes associated with our antenna, namely, 
the aperture plane and the focal plane. In principle, we can obtain one distribution from the other 
through the use of a two-dimensional Fourier transformation. This is based on the Fresnel 
approximation and will be phrased more precisely later on. First, however, we have to deal with 
some preliminaries concerning two-dimensional Fourier transforms in general. 

We start with a function of time g ( t) and its Fourier transform g ( f ) given by 

BO 

gif) = f git) e~ i2Kft dt (3.1) 

— oo 

Notice that the obvious requirement that the argument of the exponential function be dimensionless is 
satisfied by the fact that the dimension of frequency is 1/time. 

Now let us change the scenario by assuming that the function g represents a distribution of 
some entity along a line in space so that git) is replaced now by a function of x where x 
specifies position along that line. When we want to determine the Fourier transform of this 
distribution, we are faced with the choice of the substitution to be adopted for the frequency variable 
of (3.1). Goodman [2] (and many others) chooses what is referred to as "spatial frequency" measured 
in units of 1 /length. While this satisfies the above-mentioned dimensionality constraint, it is 
conceptually very unsatisfactory because usually both g and g are functions of position (i.e., g is a 
function of the position variable x, while g is a function of the position variable x) . 

We adopt here a different formulation: We preserve the dimensionality constraint by 
introducing a reference length (L) and adopting the dimensionless entity (length / L) to be the 
argument of both g and g. Specifically, we let 


X = x/L 


(3.2) 


X = x/L 

and the Fourier transform is now expressed as 

giX) = f giX) e - i2KX * dX 


(3.3) 


(3.4) 
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Generalizing to the two-dimensional case, we obtain 


g(R) = jj g(R ) e ~ i 


-i 2 n R 


R ds 


(3.5) 


where dS is a (normalized) area elment in the input plane spanned byi? . In words, the spatial 

A 

distribution ) is the two-dimensional Fourier transform of the spatial distribution g(R ) where 


R = r/L 


(3.6) 


R = r/L 


(3.7) 


It should be pointed out that the reference length L is not an arbitrary parameter. The physics 
(or other underlying theory) which yields the Fourier transformation property associating the 
functions grand g^will also yield a specific value for the reference length L. The thin-lens image 
transformation which concerns us here is a case in point. From Goodman’s formulation of this 
transformation (equation (5-15) of [2]), we infer that for an antenna of focal length F, operating at 
the wavelength X , the reference length L is simply the geometric mean of these two lengths, that is, 

L = \[xF r (3.8) 


We turn now to the precise formulation of the thin-lens image transformation which we plan 
to apply to our DO A problem. Let 


A A 

r = position vector in the aperture plane ( r = 0 on optical axis) (3.9) 

r = position vector in the focal plane (r = 0 on optical axis) (3.10) 

u = a selected field component in the aperture plane (3.11) 

u = the corresponding field component in the focal plane (3.12) 

A 


Then, the focal-plane field distribution u(R) is obtainable from the aperture-plane distribution u(R) as 
follows: 
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(3.13) 


u(R) e ~ inRi 



A 

q -i2nR'R 


dS 


where a is a constant phase. This formulation is an exact equivalent of equation (5-15) of [2]. 
Comparing (3.13) to (3.5), we note that, apart from the trivial constant phase shift a, 

A 

[u(R) e -•**«*] is the Fourier transform of u(R ) . Equivalently, we may say that u(R ) is 

A 

approximated by the Fourier transform of u(R) with the factor e 1 ^*****) accounting for the deviation 
from an exact Fourier transform. 

In our application, we start with a focal-plane distribution and look for the aperture-plane 
distribution that caused it. 

Hence we invert (3.13), getting 


u(R) = e 


-la 


//[“<*> 


, “i ICR 


1 


s i 2 it J? ‘R 


dS 


(3.14) 


This formulation provides the basis for our algorithm. In its actual implementation, however, we 
have to replace the continuous Fourier transform by a (two-dimensional) Discrete Fourier Transform 
(DFT). We proceed now to set the stage for this modification. Our starting point is the focal-plane 
array of closely packed circular horns shown in Fig. 2. Let Rj be the normalized position vector of 
the center of the aperture of horn j and let v(Rj) be its (complex) voltage output. We consider the 
selected field entity at point Rj to be proportional to v(Rj ) . 

■ In Fig. 2, we show three concentric "rings" of horns surrounding the central horn at the 

origin. Though our initial actual implementation involves only one ring, we formulate our analysis in 
terms of an arbitrary number of rings (IV) . In this general case, the total number of horns in the 
array (J) is given by 

J= 1 + 3 N(N + 1) (3.15) 


Thus, we have a set of J focal-plane voltages v(Rj) i which we want to process in order to get a set 
of aperture-plane "voltages" v(R t <) . In this context (and ignoring irrelevant scale factors), (3.14) 
transforms to 
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(3.16) 


A 

V 





A 

Q i2nR k -R } 


Equivalently, 


S'***'- *»''(*,) 

j-1 


(3.17) 


As in the case of the one-dimensional DFT, there are two important aspects of the transition 
from the continuous (3.14) to the discrete (3.16). First, the fact that we use discrete samples of the 
input function rather than the function itself means that the computed transform is periodic and could 
(potentially) be corrupted by aliasing. However, whereas in the one-dimensional case we specify this 
periodicity by a single parameter, say, the period or the frequency, in the two-dimensional case we 
have to be concerned with the size and shape of the repeating cell. 

The second aspect derives from the finite summation in (3.16), that is, the fact that we have 
only a finite number of sampling points (horns) in the focal plane. Usually, this is handled by the 
introduction of a window function. However, as shown in Appendix B, for our adopted parameters, 
the terms excluded from (3.16) are very small and may be simply ignored. Thus, the only issue to be 
addressed in the implementation of the discrete (3.16) is the effect of the aperture-plane periodicity 
imposed by it. 

In preparation for this task, we introduce now some required notational conventions. The 
focal-plane distribution of horns suggests that the nonorthogonal g ± frame shown in Fig. 2 is more 
suitable than the more conventional orthonormal h i frame shown there. If we denote the horn 
diameter by d, then any R • is expressible in terms of components in the g i frame as follows: 


Rj - j 3 — ( d/L ) (Jj^i + j' 2 ^ 2 ) (3.18) 

Note that j I , j 2 are integers ranging over the interval ( -N, N) (N is the number of rings), 
though not all such positions are occupied. For example, while there is a horn at position i? 31 there 
is no horn at position i? 3/ (see Fig. 2). 

As indicated in (3.18), we find it convenient to designate one and the same position vector as 
either a singly indexed entity or a doubly indexed entity. This, of course, implies a unique 
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correspondence between j and the pair (j 1 , j 2 ) . Indeed, in the software implementation of the 
algorithms developed here, we specifically establish such a relationship. However, in the theoretical 
manipulations in this report, it is enough to know that such a relationship exists without having to 
specify its particular form. 

It is well known that manipulations involving nonorthogonal frames (such as { g L , x ) are 
appreciably simplified by the introduction of a reciprocal (dual) set of base vectors {g 1 } which satisfy 
the cross-orthonormality condition 

9i * = 5ij U, j = 1, 2) (3.19) 


Both sets of base vectors for our case are shown in Fig. 3. Note that while any vector may be 
represented in either frame, a judicious selection of frames will yield simpler expressions. In 

A 

particular, we choose to represent an arbitrary point J? in the aperture plane in terms of the g 1 frame 
as follows: 


R = (L/d) (pj g 1 + p 2 g 2 ) (3.20) 


where p x/ p 2 are arbitrary scalars. 

We exploit this representation now to establish the aperture-plane periodicity mentioned 

A A 

above. From (3.16) we see that the dependence of v(R) on R is effected through the dot 

A 

products ’Rj , which we evaluate using the representations (3.18), (3.20) and applying (3.19) 

A 

R'Rj = PiJi + P 2 i 2 (3.21) 


A 

Now we look for a (A R) such that 


v(R + A R) = v(R) 


(3.22) 


From (3.21), (3.20), and (3.16), we infer that 

A 

A R = k 1 p 1 + k 2 p 2 


where k 1 , k 2 are arbitrary integers and 

p i = (L/c?) g 1 (i = 1,2) 


(3.23) 


(3.24) 
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Fig. 3. The reciprocal (dual) vector sets. 
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Recalling now (see Fig. 3) that 

g 1 = 2A/3 (i = 1.2 ) (3.25) 

we get for the length ofp i 

p i = (2//3)(L/d) (i = 1,2) (3.26) 

This establishes the basic periodic structure which is shown in Fig. 4. We see here that the 
aperture plane is tiled with 60° rhombi with sides of (normalized) length(2 / /3) (L/ d). The 
established periodicity means that, using (3.16) to compute v, we will find that the distribution 
of vover a rhombic tile is the same for all tiles. 

A 

Now we use this information to guide us in the selection of the set of points R k where the 
aperture-plane voltages are to be computed. Obviously, all of the J? ^ points should be constrained to 
a single rhombic tile. Furthermore, the approximation error of the Fresnel approximation (equation 
(4-10) of [2]), which forms the basis of (3.13), increases with the distance from the optical axis. 
Hence, the selected rhombic tile should be the one centered at the origin. 

Recall now that, in the focal plane, we have placed (2N + 1) points equally spaced along 
each g i axis (Fig. 2). Generalizing from common practice in the application of the one-dimensional 
DFT, one could argue that a reasonable strategy is to select the R k aperture-plane points to fully 
cover the origin rhombus with (2 N + 1) points equally spaced along each g 1 axis as shown 
(for .AT = 3) in Fig. 5, for a total number of points (J<0 given by 

K = (2 N + l) 2 (3.27) 

Note that, unlike Fig. 2 where the circles represent the circular apertures of the conical collecting 
horns, the circles here are just convenient representation artifacts to highlight the fixed spacing of the 

A 

computed aperture-plane points: there is an R k point at the center of each circle. 

Our main concern in the selection of the aperture-plane points is the avoidance of the 
detrimental effects of aliasing. As we shall see shortly, in some practical cases, this dictates a higher 
density of aperture-plane points. In view of that, we deviate from the configuration of Fig. 5 by the 
introduction of the extra scaling parameter s . Specifically, we propose to pack the K points inside a 
smaller, central, 60° rhombus given by 

(the normalized side of the R k bounding rhombus) = s j (3.28) 

(see Fig. 6) and impose the condition 

ssl < 3 - 29 > 
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Fig. 4. The tiled aperture plane. 
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Fig. 6. The adopted configuration of aperture-plane data points. 
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The importance of the scaling parameter s becomes obvious when we superimpose the 
circular boundary of the antenna aperture on the tiled aperture plane of Fig. 4. If this circle lies 
totally within the origin rhombus, then the periodicity of the computed transform becomes a purely 
academic issue and does not affect us. However, some of the K computed points will fall outside the 
circle in this case and would thus be wasted. Hence, we are better off shrinking the bounding 
rhombus, that is, adopting s < 1 in (3.28). 

If the circle extends beyond the origin rhombus, then the periodicity does manifest itself and, 
consequently, we have to deal with the effects of aliasing. It turns out that the parameters adopted in 
the simulations discussed in Section IX fall in this category and Fig. 7 shows the periodic structure 
across the full aperture for this case. 

We obviously have here a severe case of aliasing, that is, we may get significant differences 
between the "true" value of the transform and its computed value. Consider, for example, point A (of 
Fig. 7). Due to aliasing, the computed transform at this point is the sum of its true transform plus the 
true transforms at points B„ B 2 , C„ C 2 , D, and D 2 . We claim, however that, due to the symmetries 
associated with these points, the computed value still yields the true phase at point A. To see this, 

A 

denote by i? A the (normalized) position vector of point A. Apart from an irrelevant real multiplier, the 
true voltage at point A is given by 

v A = fc ) (3.30) 

Now, if we contaminate this value by adding to it the true contributions from points B, and B 2 

A 

(with position vectors R A ± p 1 ) , we get 

v AB = [1+2 cos (Lk • p 1 )] (3.31) 

and we see that the phase has not been affected. 

The same mechanism operates with the other two pairs of points yielding the more general 

result 


V ABCD 


_ A 

e J ( Lk ' x *~ vt ) {1+2 cos (Lk * p 1 ) + 2 cos (Lk • p 2 ) 


+ 2 cos [Lk * (p 1 - p 2 )] } 


(3.32) 


and our claim for point A is thus validated. 

However, when we consider point A (see Fig. 8), we find that the symmetry argument is 
applicable only to points C, and C 2 . All the other marked points are unpaired and will 
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Fig. 7. Example of a point free from aliasing. 
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thus modify the computed phase at point A. More generally, we find that all the shaded tiles on the 
boundary of the aperture in Fig. 9 represent unpaired regions of the aperture and thus cause errors in 
the corresponding areas of the origin rhombus. When we plot these "forbidden" regions of the origin 
rhombus, we find that its only region free of aliasing phase errors is the small dotted central area 
bounded by circular arcs. Trivial computations based on Fig. 9 show that in this particular case we 
should set s < 0. 156 to avoid aliasing errors. Note that not all of the unpaired regions have been 
shaded in Fig. 9. This does not affect our result because the corresponding "forbidden" regions are 
already accounted for by the regions which have been shaded. 

A 

Finally, denoting by dthe actual (unnormalized) separation of adjacent points, we have (see 

Fig. 6): 



(3.33) 


This yields the following set of aperture-plane points: 

R k = R kt .k t = W/L) [kjlff'/ff 1 ) + kj,g 2 /g 2 )] 


(3.34) 


where k t , k 2 are integers ranging over the interval ( -N, N) . Unlike the focal-plane situation 
(3.18), all such k 1 values correspond to actually computed aperture-plane points. 
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Fig. 9. Example of an aliasing-free region. 


IV. LEAST-SQUARES DETERMINATION OF A WAVEFRONT 


In Section II we found thatrz , the unit vector pointing to the source, is expressible in terms 
of the gradient of the phase distribution in the aperture plane (see (2.15)). In the absence of thermal 
noise and other types of error, three points of the distribution would suffice to determine h . To 
partially overcome the effect of noise, we provide more than this minimum number of points and are 
thus faced with a least-squares determination of h . However, instead of analyzing a least-squares 
estimation of a gradient, we prefer to deal with an equivalent— but geometrically more 
meaningful — problem, namely, the least-squares estimation of the orientation (and location) of a plane 
from a set of distances to it. 

The connection between these two approaches is provided by the following rephrasing of 

(2.7): 


n = -Vw 


(4.1) 


where 



(4.2) 


We recognize was a distance measured in the field of a plane wave along the direction of its 
propagation. Actually, this is the distance from the point with the measured phase a to the (first) 
wavefront corresponding to a = 0. 

The situation is described in Fig. 10 where we show a wavefront (orthogonal to the plane of 
the figure) at a distance a from an observer located at the origin (point O) and moving towards him. 

A 

We also show one of a set of iC points ( P k ) with position vector r k and distance w^from the 
wavefront. Obviously, 

w k = a - r k -n (1 s k s K) (4.3) 


As we shall see shortly, this time it proves convenient to normalize all lengths to d— the 
aperture-plane separation of the computed points. Hence, we replace (4.3) by 


w k /d = ( a/d ) - (r k /d) * n 


(1 <; k <; K) (4.4) 


25 



Now we infer from (3.34) that 


r k /d = = (y/3/2) fag 1 + k 2 g 2 ) 


(4.5) 


Substituting (2.13) and (4.5) in (4.4) and applying the g 1 expansions shown in Fig. 3, we obtain 


[( 1 ) fa~ k i)\ n i ~ [(>/3/2 )(k 2 + k 1 )]n 2 + (a/d) = (w k / d ) (1 s k z K) (4.6) 


A 

We have here a set of K equations to solve for the three unknowns n 1 , n 2 and (a/d) . If 
the "measured" data ( w k ) are totally error-free, then (at least in principle) any three of the above K 
equations could be used to solve for the three unknowns. If, however, w k is replaced by its 
noise-contaminated version (w k ) , then the absolute equality (=) in (4.6) is replaced by (=), that is, 
"equal in the least-squares sense" [3], and (4.6) is transformed into 

[(!) ( k 2 ~ k i)]^i ~ [(\^/ 2 )(k 2 + k 1 )]n 2 + (a/d) = (w k /d ) (1 <, k <. K) (4.7) 

where n x , h 2 and a are the least-squares estimates of n 2 , n 2 and a , respectively. 

Note that as the number of focal-plane rings ( N) increases from 1 to 3, Granges over the 
values 9, 25, and 49 (see (3.27)). Thus, we can expect significant noise cancellation in applying 
least-squares methods to the solution of (4.7). The specific method we choose is the one based on 
performing a Singular- Value Decomposition (SVD) of the linear operator involved in (4.7). We will 
show that this approach yields a very simple closed-form solution requiring neither matrix inversion 
nor any special SVD subroutines. This follows from the very special structure of the linear operator 
we are dealing with here. In the next section, we investigate this structure and apply it to the solution. 
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V. THE STRUCTURE OF THE LINEAR OPERATOR INVOLVED IN THE WAVEFRONT 

DETERMINATION 

We start with a brief overview of the SVD approach to the least-squares problem. Our 
formulation is quite general with the exception of the dimensionalities and rank, which are those 
specific to our special case. For a more detailed presentation of this approach see, for example, 
sections VII and IX of [4]. 

Assume that we are given the matrix equation 


Ag = p 


( 5 . 1 ) 


where A is a rectangular K x 3 matrix, and q, pare column matrices of three and -RTrows, 
respectively. Given A and a noise-contaminated version of p (denoted p) , we want to find the 

A. 

(least-squares) best estimate of g (denoted q) . 

The dyadic equivalent of (5.1) is 


A-q = p 


( 5 . 2 ) 


where q is a vector in an abstract three-dimensional space, p is a vector in a quite different 
K -dimensional abstract space, and A is the dyadic (linear operator) which operates on a vector in the 
three-dimensional input space to yield a vector in the X-dimensional output space. 

The association between a dyadic equation such as (5.2) and the corresponding matrix 
equation (5.1) always depends on the adoption of two specific orthonormal frames: 
one which spans the input (three-dimensional) space and one({©i}i-i) which spans the 
output ( K -dimensional) space. More specifically, the column matrix q (with elements q^) consists of 
the components of the vector q in the fj frame, 


q 5 = q-f j 


( 5 . 3 ) 


the column matrix p (with elements p i ) consists of the components of the vector p in the e i frame, 


Pi = p-h 


( 5 . 4 ) 
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and the rectangular matrix A is the representation of A in the above two frames, that is (see 
Appendix A), 


A 


ij 


* A 'f . 


( 5 . 5 ) 


Similarly, q is the vector whose representation in the({f j)j.xj frame is the column 
matrix q while p is the vector whose representation in the ({ e f rame ‘ s column matrix p. 
The fundamental SVD theorem tells us that our dyadic A is representable as follows [4]: 

A = £ 3 i U i^i (5 ' 6) 

i*i 


where the u/s are i^-dimensional orthonormal vectors, the v/s are three-dimensional orthonormal 
vectors and the s/s are positive scalars— the singular values of A (ours is a full-rank case). Dotting 
both sides of (5.6) with Vj (on the right), we see that it is equivalent to 

A-Vj = Sj Uj ( 5 . 7 ) 


Note that one can always select an arbitrary orthonormal set in the input space, operate on it 
with A and get the corresponding output set. In general, however, this output set will not be 
orthogonal. The special property of the v i set (which makes its determination nontrivial) is that, in 
addition to being orthogonal itself, the output set it generates in the output space is also orthogonal as 
implied by (5.7). Incidentally, this provides a simple test for an orthonormal set claimed to be 
the v i set. 

It is a well-known fact that, once the SVD of a linear operator (such as (5.6)) is known, its 
generalized inverse is immediately available and the solution of the associated least-squares problem 
becomes trivial (see, for example, sections VII and IX of [4]). In our case, the solution takes the 
following form: 


* E 


s/‘ VjUj-p 


i-l 


( 5 . 8 ) 
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Thus, the main computational effort in finding q is the decomposition (5.6) rather than the dot 
products involved in (5.8). 

Now we claim that, for the linear operator of (4.7), 

“ *i 

so that no decomposition is required. In order to verify the validity of this claim, it is convenient to 
introduce the vectors 


(HJS3) 


(5.9) 


z i 


A-f, 


(H I i 3) 


(5.10) 


In view of the preceding discussion, our claim (5.9) is valid if the z 1 vectors (1 £ i £ 3) are 
mutually orthogonal. Note in this context that (see (A -28) and (A-29)) 


e k‘Zi “ e k' A ' f i = A ki 


(5.11) 


so that column i of A is the representation of z i in the e k frame. Therefore, an equivalent formulation 
of the validity condition is that the columns of A be mutually orthogonal. 

To facilitate the verification of this condition, let us express the variables of (4.7) in terms of 
the new variables introduced in this section. 

n 3 = g-fj = q^ (j - 1.2) (5.12) 


a/d = <7-f 3 = q 3 


(5.13) 


w k /d = p-e. 


(5.14) 


Aic,i _ (l/2) (k 2 ~ k 3 ) 


(5.15) 
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^k,2 (v^/2) ik 2 + k j) 


(5.16) 


e *‘*3 = At. 3 = 1 


(5.17) 


We proceed now to check the three dot products in the {z d } j. x set 


N 


Z 1 'Z 3 = (1/2) J3 (ic 2 - = 0 

AT 


(5.18) 


N 

z 2 -z 3 = -(v/3/2) 53 + -fcj) = 0 

/c a , k 2 =~N 


(5.19) 




Zi'Zj = -(v/3/4) 53 (*| - *!) = 0 


k lf k 2 m -N 


(5.20) 


We conclude that (5.9) is valid and proceed to determine the remaining SVD parameters 
required for the solution (5.8). First, we infer from (5.7) that a corollary of (5.9) is 


= s i u i 


and, consequently. 

Thus, we can rewrite the solution (5.8) as follows: 


s i = 


(5.21) 


(5.22) 


£ - 53 z?fiZi'P 


(5.23) 


2=1 
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But, in view of (5.12) and (5.13), this means 



(i = 1,2) 


(5.24) 


a/d = 


z 3 -p 


(5.25) 


Of the three squared magnitudes appearing here, z 3 2 is obtained trivially as 

s 3 2 = z 3 2 = K (5.26) 

In computing the other two, we use the following well-known summation formula (Series No. (19) 
of [5]) 

J2 m 2 = (1/6 )M(M + 1) (27f + 1) (5.27) 

771=1 


which yields (using (3.15)) 

s 2 = z 2 = (1/4) £ (7c/ + 7c/ - 2k 1 k 2 ) 

k lt k^-N 


= (1/6) (2N + l) 2 N(N + 1) = (1/18) K(J - 1) 


(5.28) 


S 2 = z 2 = ( 3/4) £ (V + V + 2k i k 2 ] 

k 1 ,k 2 --N 


= (1/2) (27V + l) 2 N(N + 1) = (1/6) K(J - 1) (5.29) 

It should be pointed out that we have deviated here from the standard SVD convention 
regarding the indexing of the singular values, namely, s i s i+1 . Our hand is forced here because 
the ranking of the s/s (by value) varies with N. 
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(5.30) 


Applying the values obtained here, we get the very simple final results: 

A N 

n x = {9/ [ K(J - 1) d] } £ {k 2 - k x ) w kiikj 

k^kj—N 

n 2 = -{3j3/[K(J - l)d] } X) {k 2 + (5 - 31) 

a.(l/Jof 'A < 5 - 32 > 

Jc-1 

Note that the equations for n 1 , h 2 may be expressed in terms of templates applied to 
the w k kj distances, that is, each w ki ki is multiplied by the template value corresponding to its 
position in the aperture plane and the estimated parameter is the sum of these template-weighted 
distances. These templates are shown in Figs. 11 and 12, which display very distinct simple patterns. 

These patterns suggest that, as long as we restrict the R k points to a small fraction of the 
basic-period cell, placing the K aperture-plane points inside a square (rather than a rhombus) is 
probably an equally good alternative to the rhombic arrangement. 

Though we now have a simple algorithm which yields the best 0 east-squares) estimate of the 
DOA, its practical significance is rather minimal unless we complement it with a second algorithm 
which will provide us with information on the uncertainties associated with this estimate. It turns out 
that this is a much more challenging undertaking; it will occupy us in the next three sections. 

In preparation for this task, we point out two results which, though of no direct interest, are 
required in the development of the second algorithm. The first one is the (trivial) value for the 
estimate of a —the distance to the reference wavefront (5.32). The second one is the complete 

A 

vectorial representation of the least-squares estimate of n , which we denote n . Obviously, 

n = n 1 h 1 + n 2 h 2 + 

where {h is the antenna-fixed frame (see Fig. 1). 
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VI. CONFIDENCE BOUNDS FOR THE DOA DETERMINATION (PRELIMINARY 
CONSIDERATIONS) 


The DOA may be represented as a point (D) on the unit sphere (the intersection of n with the 

A A 

unit sphere). Having solved for n , we are now in possession of the location of point D, which is the 
best (least-squares) estimate of point D. Now we want to know the reliability of this estimate. 

, a 

Specifically, we want to know the size (and shape) of the region around D where the true point D is 
likely to be found with a prescribed probability (say, 99%). To attain that goal, we have to start with 

A 

the statistical properties of the least-squares vector from which n is derived, namely, the least-squares 

A 

solution vector q . 

First, we fecall that the very nature of the least-squares processing implies 

A 

E{q) = q ( 6 . 1 ) 


where Sis the expectation operator. Expressing this in terms of the first-order central moment 

A 

of q, we have 


- g) = o 


( 6 . 2 ) 


r 


% 


A 

To attain our goal, we have to evaluate the second-order central moment of q , namely. 


Y = E{ (g - g) (g - g) } (6.3) 

Note that Y is the covariance dyadic whose relationship to the (more familiar) covariance 
matrix Y (with elements is expressible as follows: 

y t Y a i i ~ f i < 6 - 4 > 

i,J m i 


To evaluate Y , we have to complement the expression for g (5.8) with a similar expression 
forg, namely. 


3 

i 

i-1 


g = E s i _1 ^i u i‘P 


(6.5) 
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This is simply the special case of (5.8) in which there is no noise corrupting the data vector p . The 
only validity condition (for both equations) is the requirement that A should have full rank (3) or, 
equivalently, that all three singular values be nonzero. This is, indeed, the case as evidenced by 
(5.26), (5.28), and (5.29). 

Subtracting (6.5) from (5.8) and recalling (5.9), we obtain 


A 


q - q = 


3 



S i _1 f i “ P) 


( 6 . 6 ) 


The vector (p - p) appearing here represents the aperture-plane noise, which is the cause of the 
uncertainty in q . Let us denote 


< 6 - 7 > 

In terms of components in the e i frame, this means (see (5.14)) 

Ci = ©i-C = ("i - Wi)/d ( 6 -8) 

Now we apply (6.6) and (6.7) to the evaluation of Y (6.3), getting 


- E 

1 . j-i 


-l 


f i u i 


•£{((} 


u J f J 


(6.9) 


In view of (6.4) and (5.21), this means 


= (6.io) 

lj Zi z i 

The only unknown entity here is E{ C C } , which is the covariance dyadic of the noise contaminating 
the data vector of the least-squares problem. In the usual least-squares application, one assumes the 
corresponding covariance matrix to be diagonal and, based on that, simplifies the expression. In 
our case, the situation is quite different because the noise is injected into the system not in the 
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aperture plane where the least-squares algorithm is applied, but in the focal plane. Thus, while the 
focal-plane noise is, indeed, uncorrelated, the processing that generates the aperture-plane distribution 
leads to highly correlated aperture-plane noise. Our immediate goal, then, is the evaluation of the 
covariance dyadic of this noise, which we take up in the next section. 


VII. THE COVARIANCE OF THE APERTURE-PLANE NOISE 


As outlined in the last section, our task here is to evaluate the dyadic E{ ( ( } or, 
equivalently, the set of scalars E{ C k C m } for all k, m . The parameters which determine C * are 
illustrated in Fig. 13, which shows the k -th noise-free complex voltage with magnitude a k and phase 

A 

angle a* . This voltage is contaminated by the added complex noise voltage N k (with 
magnitude N k ) , which leads to the phase error as shown. From (6.8) and (4.2), we infer 

C k = (r\ k /2n) (X/d) (7.1) 


so that 

WjcC*} = [X/(2nd)) 2 E{l) k n\ m } (7.2) 


and we concentrate our attention initially on £'{Tij t r| ni } . 

Though we are dealing here with complex scalars, it is convenient to regard them as 
three-dimensional vectors (over the field of reals) in a space with an orthonormal basis {CjJj.j, as 
indicated in Fig. 13. Note that, in this isomorphism, the vector c 1 is isomorphic to the real number 1, 
the vector c 2 is isomorphic to the imaginary number i, and the vector N k is isomorphic to the 
complex number N k . The motivation for establishing this isomorphism with three-dimensional 
vectors rather than two-dimensional vectors will become clear shortly. 

The computation of £{ri A .ri m } in the general case presents serious mathematical difficulties. 
However, if we constrain our analysis to the practically significant case of large SNR, the problem 
becomes manageable. We propose, then, to pursue this initial analysis under the constraint 

N k << a k (7.3) 

This means that, even in the worst-case orientation, we have 

ri* = tan q*. (7.4) 


so that 


n* * 


N k sin ijr*. _ N k sin 


a k + N k cos ty k 


(7.5) 
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Fig. 13. The aperture-plane phase-error geometry. 
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But this is expressible as a vector product, namely. 


Hjc^3 


a k x N k 


(7.6) 


Hence, we can express the product (T^ru) as a dot product: 





(7.7) 


Applying now elementary vector manipulations, we reduce this to 


(^A) 2 (n*n„) - (** -a m )(N k -N m ) - a m ‘N k N m ’a k (7.8) 


and finally. 


( 3 k a ,) 2 - ^£{^4 •** (7 . 9 ) 


To evaluate the expectations appearing in (7.9), we have to trace the evolution of N k from the 
noise generated in the focal plane. The motivation here is that we can safely make reasonable 
assumptions concerning the statistical properties of the noise added to the focal-plane horns’ outputs. 
Let the complex noise added to the output of the j - th horn be rr, . Applying (3.17), we find that 



e inSj-i2& k - 


(7.10) 


Applying now the isomorphism introduced earlier, we may represent the complex noise rij by the real 
vector rij and (7.10) will have to be expressed in terms of rotation operators. Specifically, the 
complex scalar is isomorphic to the real vector T’(c 3 , P) ~n j where 


| 
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T(p , P) = a rotation dyadic which, on operating to the right, rotates a 

vector by the angle P about an axis given by the unit vector p (7.11) 

Though there is a relatively simple explicit formulation of T( p , P) in terms of its p , P parameters, 
we do not need this here. For our derivations, it is sufficient to realize that T(p , P) is a 
length-preserving linear operator (represented by an orthogonal matrix) and thus satisfies 


T' 1 (p-,p) = T(p, -p) = T(p , P) 


(7.12) 


Applyingjhese ideas to (7.10), we get 




T[c z ,nR j '(2R k - Rj)]-n. 


(7.13) 


To simplify the notation, we denote 


A 

T kj s T[C3'*Rj‘( 2R * ~ R j)] 


(7.14) 


so that 


u vs 

~ ^ "Tic j' n j ~ ^ ^ j'^k j 
j - i J-i 


(7.15) 


Now we apply (7.15) to evaluate the two expectations appearing on the right of (7.9). 
Specifically, 


J S _ 

E{N k -N m ) = 52 (7.16) 

j.s - 1 
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(7.17) 


J _ £ 

= E 

j , S*1 


To proceed further, we have to utilize now the statistical properties ofn y, the noise component of 
the output of the j-th horn. We adopt the following two assumptions concerning these noise 
components: 

(1) All noise components (in-phase and quadrature) of all the horns’ outputs are totally 
uncorrelated. 

(2) The in-phase and quadrature noise components of the j-th horn have the same 
power (of) . 

To facilitate the mathematical formulation of these assumptions, we express n j in terms of its 
components as follows: 


n j = E n ju c u (7-18) 

u-i 

Our two assumptions are equivalent to 

E{n ju n sv ) = o 2 b js b uv (7.19) 

We apply this equation to (7.17) first, starting with the evaluation of 

= 

where I is the identity dyadic in the space spanned by the frame {Cj} 3 iml . Using this result in 
(7.17), we obtain 


E E { n ju n sv) c U C , 


U, V=1 




2 

E 

U = 1 

- ^ 3 ^ 3 ) 


(7.20) 
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W k x„) - E °/ - c 3 c 3 yT mj 

i-i 


j- _ ^ \ 

Vj=i 


- C 3 C 3 


X>/ 


( 7 . 21 ) 


j=i 


We have used here the obvious fact that a rotation operator applied to the rotation axis (c 3 ) leaves it 
unchanged. Resorting now to the T k j definition (7.14), we can replace the product of the two 
rotation dyadics we have here by a single rotation dyadic as follows: 


T kj’ T mj = T[C i ,-%Ry(2R ia - Rj) + uRj "(2R k - R j)] 

— A A 
= , 2u(R k - 


( 7 . 22 ) 


The angular argument appearing here turns out to be an important system parameter. We denote 
it Y *,0, that is, 

y kmj ~ ^ - R^’Rj (7 


so that 


T]cj ‘ T m j T’(C 3 , Y kmj) 


( 7 . 24 ) 


and, finally, 


E {*W 


j 


E °/^( C 3'Yicmj) 


^3^3 E ° 


j-l 


2 


( 7 . 25 ) 
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We turn now to the other expectation, namely (7.16). Applying (7.19) to it, we obtain 


E {N k 'N m } - ^2 E E { n ju n s-^ C u'^kj ' ^ns C v 

j,s=l u,v-l 

J-l u-l 


(7.26) 


Repeating the manipulations involved in (7.22), we get in the present case 


T kj -T mj , -y J onj) 


(7.27) 


so that 


E {N k 'N m } = J2 E °/ [ 5 U * r ( c 3' -Yjw)*c u ] 

Jf-l ' 


U-l 


(7.28) 


Having computed the two required expectations, we go back now to (7.9) and substitute 
(7.25) and (7.28) in it. Recalling that the vectors a k , a m are orthogonal to c 3 , the result is 


J 2 _ 

(a Jt a /n ) 2 E{r 7jc r 7m } = (i^a^ E [^‘^( c 3' "W c u] 

j-1 U*1 

- E °1 K' ? ( 5 i' < 7 - 29 > 

J*1 


We have here two quadratic forms involving a rotation dyadic. Since all the vectors in these 
quadratic forms are orthogonal to the rotation axis (c 3 ) , the geometric interpretation of these 
quadratic forms is very simple. Consider the second quadratic form while reviewing Fig. 13. 

Vector a* has magnitude a. k and forms an angle a k with c 3 . Operating on it with T(c 3 , y j^) yields 
a vector of magnitude a k forming an angle (a k + Y^) with c k . Finally, dotting this vector 
with a m yields the result 
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a «*^,Yj ^ya k = a k a m cos (a* - a m + y^) (7 . 30; 

Applying the same reasoning to the other quadratic form yields 

2 

E C U 'T(C 3 , -y jon-jj'C u = 2 COS Yj fcnj- ,-j 31 - 

u=l 

These results simplify (7.29), yielding (after some trivial trigonometric manipulations) 

j j 

(ajca m )£{n*n4 = cos («*-«„) V Of cos ykmj + sin(a Jt -a fl) ) £ 

j=i j-i 

( 7 . 32 ) 

The two sums appearing here are system parameters involving geometric and noise entities, but these 
sums are independent of the signal. For a given system configuration, they are constants which have 
to be evaluated only once. We introduce now the following symbols to represent these sums: 


= E a / cos [ 2 ~ Z»)- R j] 

J m 1 


3 ** = E o/sinp*(JZ* - R m )-Rj] 

j- i 


Note that we have replace! y ^ by its value (7.23) to bring out the explicit dependence on the 
geometry. In terms of these symbols, our final result is 

= [l/(a Jt a jn )][C Jtm cos(a Jt - eg + sin (cc k - eg] (7.35) 


This formulation allows us to handle the case where different amplifier types are used for different 
horns or to accommodate variations in the noise performance of supposedly identical amplifiers. 
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One cost-effective design would use identical (cheaper, higher-noise) amplifiers for all horns 
except the central one [6] . In this case, we can take advantage of the symmetry of the distribution of 
the horns in the focal plane. Note that all horns other than the central one occur in pairs at 
locations ±R j (see Fig. 2), and thus their contributions to cancel out. The central horn makes no 

_ a 

contribution either since the corresponding R j is zero. Thus, vanishes in this case and we get the 
simpler result: 


spun*} = 


A 



COS (OCfc - eg 

® k^m 


(7.36) 


Finally, if all the amplifiers are identical with the common noise power o 2 , then 

4» = ° 2 cv 

where C ^ is a purely geometric factor given by 


'km 


(7.37) 


A A 


Cjbn ■ E cos [2n(R k - R m )-Rj] 
7-i 


(7.38) 


and the result simplifies to 


s {n*n m } = 


(7.39) 


where 


_ ^ cos (a k - eg 
— 


ton ^ton 


l 


(7.40) 


This is the formulation we adopt for the Yy expression (6.10). Hence, applying (7.2), we get 


S {C*C4 = [oA/(2ud)] 2 Hj bn 


(7.41) 


and the (previously) unknown factor of (6.10) is now expressible in terms of the ’s as follows: 


WO = [oX/(2nd)) 2 £ 


(7.42) 


where if is the dimension of the output e -space. Equivalently, 
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(7.43) 


i?{CO = [o\/(2nd)?» 


where 


K 

1 ^ “ V'km&k&m 

k,m*l 


(7.44) 


Equation (7.42) still contains one unknown, namely, o 2 . Let us elaborate: In principle, one 
could determine a 2 by direct measurements. However, in practice, this raises a lot of problems such 
as calibration, variation of o 2 with time and environmental conditions etc. Instead of that, one could 
obtain an estimate of o 2 as an additional output of the least-squares processing. Specifically, with the 
proper processing, the vector of residuals (v ) yields an estimate of o 2 . Let us examine this 
processing. 

The residuals vector (v) is defined as follows: 


A — A 

v ■ p - A -g 


(7.45) 


Substituting for A and gfrom (5.6) and (5.8), respectively, we get 


V = £ u.Ui-p 


1*4 


(7.46) 


But, for i > 3, 


u i 'p = 0 


(i > 3) 


(7.47) 


Hence (6.7), 


A 

Ui’P s U i< 


(i > 3) 


(7.48) 
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and (7.46) reduces to 


K 





( 7 . 49 ) 


Using this form of v , we find that E{v 2 } is expressible in terms of o z . Specifically, 


E{v 2 } = V Ui-£{ CO *Ui 

1^4 

= [aA/(2ird)] 2 g Uj-Ji-Ui 


( 7 . 50 ) 


In principle, (7.50) yields a 2 . In practice, we would like to avoid the computation of the 
(nonunique) set { u J. 4 This can be effected as follows: 


K 

JC ) 


3 ) 

E = 

E 

- 

E “i*l** U i 

1-4 



(1=1 / 


( 7 . 51 ) 


Note that the first sum on the right is the trace of |i ( tr { p) ) . But the trace of a linear operator is 
an invariant which has the same value in any orthonormal frame. In particular, 


tr (ji) = g e^ix-e, * 


i-l 


This removes the last obstacle to the computation of o 2 and yields: 


[ok/ (2 ltd) ] 


£{v 2 } 


K ^ 

E 1 *ii 


E 


( 7 . 52 ) 


( 7 . 53 ) 
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The covariance dyadic E{ f C } we set out to derive is thus finally given by 





( 7 . 54 ) 


where ji is given by (7.44) and (7.40). 

Note that (7.40) is formulated in terms of the noiseless parameters a k , cc k (see Fig. 13). 
Though these are not directly available to us, they are the statistical means of the corresponding 
noise-contaminated entities, which are available. One possible way of utilizing this fact in the 
application of (7.54) is as follows: We partition the input data into a large number M(say, M= 100) 
of independent data sets, apply the least-squares solution to each one of these sets, and obtain the 

A 

DOA estimate from the mean of these M individual solution vectors iq) . The covariance dyadic of 
this mean solution vector is just l/Afof the covariance dyadic Y (6.3) of a single solution vector. 
Now, in evaluating required by the single-solution Y (see (6.10)), we exploit the fact that we 

have at our disposal a large number (Af) of independent sets of noisy a k s and a k s: We simply 
compute their means (over the M sets) as good estimates of the noise-free a k s and a^’s required 
(indirectly) by (7.54) (see (7.40)). Similarly, for the factor 2?{v 2 } of (7.54), we substitute the mean 
of the M different values of v 2 obtained from the M least-squares solutions. 

As we shall see in Section IX, this approach yields reasonable agreement with results based 
on direct covariance computations. 
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VIII. CONFIDENCE BOUNDS FOR THE DOA DETERMINATION (FINAL FORMULATION) 


Having determined the noise covariance in the aperture plane, we resume now the 
development of the DOA confidence bounds. Our first step is to update the expression for the output 
covariance matrix (6.10) by substituting (7.42) in it: 


Y U 


[ok/ (2nd) ] 2 


z?z? 

z i z j 


E V** ( Z i' e k) { Z j' e m) 


k.m* 1 


( 8 . 1 ) 


All the variables appearing here are now available to us. Specifically, 

[a A./ (2nd) ] 2 in equation (7.53) 
z? in equations (5.26), (5.28), and (5.29) 

K in equation (3.27) 

in equations (7.40) and (7.38) 
z j • e k in equations (5. 15) through (5.17) 

Thus, we have here all the information required to encode (in software) the computation of the output 
covariance matrix Y. This means that the central moments of the probability density function of the 

A 

vector q (Section VI) up to (and including) the second order, are now in our possession. Nothing has 
been said, however, about the higher moments of this density function. Yet, in order to obtain the 
desired DOA confidence bounds, we have to adopt a specific reasonable density function consistent 
with the computed central moments. Our choice is the three-dimensional Gaussian density function 

A 

with the computed Fand E{q} (= q) as its parameters. Specifically, denoting 

J. < 8 - 2 > 

Z = Y~ x ( Z = Y' 1 ) (8.3) 


we adopt the probability density [7] 
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G, (x) = 


( 8 . 4 ) 


N (2n) 3 

Obviously, the surfaces of constant G 3 ( x) in the three-dimensional x-space are characterized 


by 


x-Z-x = constant 


( 8 . 5 ) 


But this is an equation of an ellipsoid centered at the origin. This leads to ellipsoidal confidence 
bounds for the vector q . However, whereas the first two components of q (n x , n 2 ) prescribe the 

A 

desired DOA, its third component ( a/d ) is of no interest to us and unnecessarily complicates the 
interpretation of the results. 

Such a situation is encountered in many other least-squares problems where it appears to be as 
unavoidable as it is in ours. The root of the problem seems to lie in the fact that, in some cases, the 
engineering parameters we try to extract are only part of a larger set of interrelated physical 
parameters. It is this interrelationship which dictates that the smaller subset of engineering parameters 
is retrievable only in the context of solving for the larger set of physical parameters. We claim, 
however, that though we cannot prevent the occurrence of such a scenario, we can cancel its 
undesirable effects. Let us examine how this is brought about. 

G 3 (x 0 ) provides the infinitesimal probability dP 3 (x 0 ) of finding x in an infinitesimal 
(three-dimensional) region in x-space around point x 0 , that is, 

dP 3 (x 0 ) = G 3 (x 0 ) dx x dx 2 dx 3 (8.6) 

However, we are interested in the probability of an (x lt x 2 ) combination regardless of the value of 
x 3 . To facilitate the formulation of such a probability, we introduce now 

(8 - 7) 

i=l 


Thus, 


A - . 

dP 3 (x 0 ) = dP 3 (x 0 + x 3 f 3 ) 
and the probability that concerns us is 


( 8 . 8 ) 
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(8.9) 


dP 2 (x 0 ) = f dP 3 (x 0 + x 3 f 3 )dx 3 


f G 3 (x 0 + x 3 f 3 )dx 3 


f dxi dx 2 


This is equivalent to the introduction of a new two-dimensional density function 


G 2 (x) 


= Jg 3 (x + x 3 f 3 ) dx 3 


( 8 . 10 ) 


A 


A 


which provides the infinitesimal probability dP 2 (x 0 ) of finding x in an infinitesimal 

A A 

(two-dimensional) region inx-space around point x 0 , that is. 


dP 2 (x 0 ) = G 2 (x 0 ) dx 3 dx 2 


( 8 . 11 ) 


With the adopted Gaussian distribution G 3 (x) , the computation of G 2 (x) is trivial: All we have to 
do is replace the three-dimensional covariance matrix Y by its upper left two-dimensional 
submatrix Y (see, for example, result 7.4.3 of [7] concerning the Gaussian marginal distribution). 
Thus, introducing the inverse 

Z = Y - 1 ( 8 - 12 > 


A 

and regarding it as the representation of dyadic Z in the f ± frames, we get 


G 2 {x) = 


N (2n) 


-i 5 


A * A 

Z-x 


(8.13) 


Having derived the explicit expression for G 2 (^) / we are now in position to compute the 

A A 

probability of finding x inside any region in x-space by simply integrating (8. 1 1) over that region. 
This computation becomes almost trivial when we select thex-space integration region to be an 
ellipse determined by the spectrum of Z and perform the integration via the corresponding elliptic 
annuli. Let us consider this in detail. 
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We start by tracing a curve in x-space along which G 2 (x) has a constant value. Specifically, 
we choose the curve which satisfies 


-X'Z'X = t 

2 


(8.14) 


where t is a fixed parameter. The shape of this curve is readily obtained from the spectrum of Z as 
follows: Note that the symmetry of Y implies that Z is also symmetric. Applying the general SVD 
representation theorem (5.6) to the symmetric case, we see that Z is representable as 

! - E (815) 

i- 1 


But, since the v i ’s are orthonormal, this means that 

A 

Z'V ± = X i v 1 


(8.16) 


and we identify the X d ’s as the eigenvalues of Z and the v i ’s as the corresponding eigenvectors. 
Substitution of the (8.15) representation in (8.14) now yields 






2t/X x 2t/X 2 


= 1 


(8.17) 


where 


h = x*v i 


(8.18) 


In other words, the selected curve is a member of a family of concentric ellipses with parameter t . 

All members of the family have their axes along the eigenvectors v 1 , v 2 while their semiaxes depend 
on the parameter t as follows: 

length of i-th semiaxis = ~t/X i (8.19) 

As we increase t, the size of the ellipse increases. Its orientation and shape, however, are not 
affected, being independent of t . 
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The area enclosed by any of these ellipses (A) is given by 


A = 2 = 2nt/f\f\ 


( 8 . 20 ) 


so that the differential area of the infinitesimal annulus along the curve is given by 


dA = ( 2 n/J ) dt 


( 8 . 21 ) 


Therefore (see (8.13) and (8.14)), the probability of finding x in this annulus is (e t dt) , and the 

✓s 

probability of finding x inside an ellipse with parameter t (denoted Pi t) ) is given by 

c 

Pit) = f e~ s ds = 1 - (8.22) 

0 


The interpretation of this result is quite straightforward. Let us illustrate it with an example: 
Assume that we have applied the least-squares solution and our best estimate of the DOA is 

A 

represented by point D on the unit sphere (see the first paragraph of Section VI). Now we want to 
find a 99% confidence region around this point. From (8.22) we find that the value of the 
corresponding t is 4.6. Substituting this value in (8.19), we find the semiaxes of the confidence 
ellipse, while its orientation is given by v x , v 2 , and its center is given by the least-squares solution. 
Strictly speaking, this ellipse is traced in the (x, y) plane of Fig. 1 and is the projection of the 
confidence region traced on the unit sphere. However since, in the usual circumstances, the center of 
this ellipse is very close to the origin and its size is very small compared to 1, this projected region is 
almost identical with the confidence region on the surface of the unit sphere. 
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IX. SOFTWARE VERIFICATION 


The main result established in the preceding sections is the development of two algorithms to 
process the outputs of the focal-plane horns. The first yields the best (least-squares) estimate of the 
DOA, while the second determines the confidence region around this estimate. Both algorithms have 
been implemented in software and are incorporated in a quite flexible program which provides easy 
selection of a large number of various operational parameters and allows the controlled injection of 
noise. 

Though this program is capable of extensive testing of all aspects of the algorithms, so far we 
have carried out only the minimal number of tests sufficient to establish their validity. Since real data 
are not yet available, we have applied this program to simulated data provided to us by P.W. Cramer 
of JPL. These data refer to an idealized version (no deformation) of the 34-m antenna at DSS 13 of 
the Goldstone Deep Space Communications Complex. The assumed frequency is 32 GHz with a 
corresponding 3-dB loss at a deviation of about 8.4 millidegrees from boresight (the spherical 
coordinates) [8], The data set consists of the outputs of conical horns of 4.4-in. (11.2-cm) diameter 
in configurations of up to three rings (see Fig. 2) and for incident waves arriving at zero to 12 
millidegrees off boresight (in steps of 2 millidegrees) with several selected values of the spherical 
coordinate <J) . 

Our first concern was the behavior of the DOA algorithm in the absence of noise. Figure 14, 
which was obtained for the three-ring array (N = 3) with<|> = 0° , is a reasonably good 
representative of the results obtained in all other cases ( N = 1,2; <|> = 90° , 180° , 27 0° ) . We 
see here a deviation of the estimated 0 from the true 0 which increases rapidly with the value 
of 0 . We shall shortly introduce a plausible explanation for this behavior. However, it should be 
stressed at the outset that this deviation does not affect the viability of the algorithm. The basis for 
this statement is simply the fact that the observed deviation follows a well-defined, smooth curve and 
could thus be easily corrected via a "calibration curve" such as that of Fig. 14, regardless of the 
details of the error mechanism. 

In searching for the possible cause of the deviation, we examined first the algorithm which 
translates the aperture-plane phase angles to the first two components oih, the unit vector pointing to 
the source. This algorithm, which is given by equations (5.30), (5.31), and formulated in terms of 
templates in Figs. 11 and 12, is so simple that it can be easily tested with a hand calculator. Because 
of the linearity of the algorithm, it is sufficient to test it with just two properly chosen inputs. The 
inputs we have selected for this purpose are 4> = 0° (n 2 = 0) and<|> = 90° ( n 1 = 0). 

According to (5.33), these correspond to 
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n 2 = sin 0 (<J) = 0°) 


( 9 . 1 ) 


n 2 = sin 0 (<J> = 9 0°) 


( 9 . 2 ) 


Trivial trigonometry now shows that all values of ( w k / d) are multiples of sin 0 . Application of 
the templates of Figs. 11 and 12 to these computed values yields the correct results ((9.1) and (9.2)) 
thus verifying the validity of the algorithm. 

We conclude that the problem of Fig. 14 must arise before the application of this algorithm, 
namely, in the transformation of the data from the focal plane to the aperture plane. In other words, 
the suspect is the Fresnel approximation and the thin-lens Fourier-transform property based on it (see 
Section m). Indeed, the very nature of the Fresnel approximation implies that it is applicable only to 
systems in which the angular deviation of rays from the optical axis (0) is small. Now, the actual 
numerical value of 0 which is considered small in this context, may differ quite a bit from system to 
system. One might argue, then, that Fig. 14 simply tells us that, for the DSS 13 antenna, the validity 
bound for the Fresnel approximation is about 5 millidegrees. The fact that the transition out of the 
Fresnel approximation region is as smooth as indicated in Fig. 14 is certainly a welcome result. As 
already pointed out, this allows us to extend the application of this approximation beyond this 
boundary. 

It should be stressed that the above explanation is just a plausibility argument. The rigorous 
mathematical analysis required to substantiate it has not been carried out. 

So far we have considered the behavior of the overall DOA algorithm with respect to the 
spherical coordinate 0 . Its behavior with respect to the spherical coordinate <J> is quite satisfactory. 
This fact lends further support to our explanation of the 0 behavior. 

We turn now to the second algorithm, namely, the one that determines the confidence region 
around the DOA estimate. Figure 15 is representative of what is obtainable with this algorithm. We 
show here a projection onto the (x, y) plane of a very small region of the unit sphere around 
the z axis. The tick marks along the axes correspond to the spherical coordinate 0 (on the unit 
sphere) measured in millidegrees. Point D represents the true position of the unit vector pointing to 

A ^ 

the source (n) . Its spherical coordinates are 0 = 0°. 004, <|> = 0° . Point D represents n —the 
estimate of Bunder the following conditions: 
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(9.3) 


(signal power) / (noise spectral density) = 40 dB-Hz 
(number of samples used in the DOA estimate) = 100 

► 

(integration time for each sample) = 0.2 sec 
(number of focal-plane rings, N) = 1 

* 

A 

The ellipse shown centered at point D is the 99% confidence region computed by the second 

A 

algorithm for this case; that is, it is the region around point D inside of which point D should be 
found with a probability of 99%. Note that point D, indeed, lies inside this ellipse. It should be 
pointed out that, though the region looks circular, it is in fact an ellipse with axes tilted with respect 
to the ( x , y ) axes. 

To check the validity of this algorithm, we carried out a direct computation of the output 
covariance matrix y(as opposed to its computation based on the theoretical (8.1)), using a much 
larger number of samples (2000). In Fig. 16, we show the confidence region corresponding to this 
direct computation superimposed on the theoretical confidence region of Fig. 15. Though the two are 
not identical, it is gratifying to note how close they are. 

The practical operational significance of a result such as that in Fig. 15 is shown in Fig. 17 
where we have added a circle of 3.15-millidegrees diameter to represent the locus of points where the 

A 

loss would be 0. 1 dB if the antenna is pointed at the estimated source direction D [8]. We see that the 
99% confidence region lies well inside the 0.1-dB circle. Hence, there is at least a 99% probability 
that the loss due to the error in estimating n is less than 0.1 dB. Of course, when more detailed 
beam-pattern measurements become available, the single 0.1-dB circle would be replaced by a family 
of concentric circles for smaller losses, and a finer estimate of the probable bounds on the pointing 
loss would be obtained. In this case, one could determine a sequence of confidence regions 
corresponding to a sequence of signal integration times, thereby getting a clear presentation of the 
trade-off between signal integration time and pointing-error loss. 

As we have already pointed out. Figs. 14 and 17 are typical of results obtained with other 
parameter combinations. Reviewing all these results, we conclude that both algorithms are viable. 
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Fig. 17. The pointing-error loss (N = 1) . 
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X. 


CONCLUDING REMARKS 


Consider the pointing of a highly directional antenna at a desired source. Assume now that the 
signal has been properly acquired so that, initially, there is no pointing error. As the source moves 
across the sky, the tracking system moves the antenna in an effort to keep the source close to the peak 
of the antenna pattern. A commonly used tracking system (CONSCAN) keeps the antenna in constant 
motion so that its axis (boresight) traces a small circle (on the unit sphere) centered at the nominal 
position of the source. When the actual position of the source drifts away from the center of the 
traced circle, the variation of the received power along the circle provides the information required to 
reposition the antenna and bring the source back to the center of the circle. 

Note that with this scheme, the antenna is always pointing a little off the source, thus 
incurring a signal loss. To minimize this loss, the diameter of the scanned circle (expressed in 
degrees along a great circle on the unit sphere) should be much smaller than the 3-dB width of the 
antenna pattern. For highly directive antennas, this may lead to such a small circle that the 
mechanical part of the antenna-pointing system is incapable of tracing it with the required precision. 

Against this background, consider now a tracking system based on the ideas developed here. 
We install in the focal plane of the antenna an array of horns (with their associated receivers) 
surrounding the customary single central collecting horn. Now, instead of constantly moving the 
antenna in a circular pattern around each point of its predicted tracking trajectory, we just process the 
horns’ outputs, thereby getting a continuous monitoring of the pointing error. We issue a command 
to deviate from the predicted trajectory only in response to a sufficiently large detected pointing error. 
Nowhere in this system do we face the constraints of the mechanical pointing system. 

Presumably, a pointing system using a focal-plane array has not been previously implemented 
because of the cost of the required additional hardware. As we have already pointed out, this 
argument does not apply to the case of our immediate concern because the hardware referred to is an 
integral part of the deformation-compensation system. However, when a very large single-dish 
antenna is considered, the investment in a focal-plane array would probably be only a small fraction 
of the total cost of the antenna. In this case, then, the antenna-pointing system proposed here would 
be a viable approach even for a nondeformed antenna. 

We have shown that the DOA algorithm developed here is equivalent to the solution of the 
purely geometrical problem of finding the best (least-squares) estimate of a plane from the noisy 
measurements of its distances from a large set of given points. In Appendix B we came across a 
different purely geometrical equivalent of the DOA problem. This alternative formulation just 
replaces the plane of the adopted formulation with a point. Specifically, we found in Appendix B that 
the DOA problem is also equivalent to finding the best (least-squares) estimate of a point from the 
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noisy measurements of its distances from a large set of given points. Though we have not followed 
through on this alternative approach, it seems worthwhile to examine the algorithm it generates and 
compare it to the one developed here. 

Another subject which merits further study is the effect of the finite areas of the focal-plane 
horns. Our analysis regards the horns as idealized point-samplers of the focal-plane field. In other 
words, we assume that the voltage output of a horn is proportional to the field at the center of its 
aperture. In reality, this voltage is proportional to the integral of the field over the aperture. Thus, 
our assumption is equivalent to approximating the field at the center of the horn aperture by the 
average of the field over the full horn aperture. For sufficiently small horns, this is a reasonable 
approximation. Though the simulation results of Section IX seem to provide an indirect validation of 
this assumption, the subject merits closer scrutiny. 

Finally, it would be very worthwhile to study the derived output-covariance expression (8.1) 
and apply the insight gained to the selection of optimal system parameters. In this context, it should 
be pointed out that the adopted set of aperture-plane points of (3.34) underlying (8.1) is just one of 
several alternative candidates. Optimization analysis would probably prove fruitful here too. 
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APPENDIX A: REAL DYADICS 


A dyad (D) is an entity represented by two vectors as follows (recall that we use overbars to 
denote vectors): 


D = ab 


(A-l) 


There is no implied operation between a and b . The meaning of the dyad is defined in terms of 
what it does to other vectors. Thus, 


D-h = (ab) -h = a(b-h) 


(A-2) 


that is, the result of (scalarly) postmultiplying D with h is the vector a multiplied by the 
scalar (b'h) . Note that, in general, the vectors a ,5 would belong to different spaces and would 
have different dimensionalities. The only constraint in (A-2) is the obvious one that A and b belong 
to the same space. 

Equation (A-2) shows D to be a (rather limited) linear operator operating on vectors in 
the 5-space to yield vectors in thea-space. This same dyad can also operate in the opposite direction. 
Thus, 


y-D = y-(ab) = (ya)5 


(A-3) 


Alternatively, this "opposite-direction" operation may be implemented using the concept of 
transposition. Specifically, the transpose of the dyad D (denoted D) is given by 


D 


= ba 


(A-4) 
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Hence, (A-3) may be expressed, equivalently, as 


D-y = ( ba ) *y = b(a-y) 


(A-5) 


Though a dyad is a very limited linear operator, a linear combination of dyads (called a 
dyadic) can represent any linear operator. Specifically, given the set of dyads {//_*} " ml , we can 
construct the dyadic H as a linear combination of its elements as follows: 


H = £ Yi^i 


(A-6) 


where the y i ’s are scalars. 

It should be pointed out that the above dyadic H may also be represented as a different linear 
combination of a different set of dyads. In other words, representation (A-6) is not unique. There is, 
however, a unique canonical representation of dyadics which is very useful and provides a powerful 
tool in manipulating linear operators in general. Specifically, an arbitrary dyadic (linear 
operator) A of rank k is always representable as follows: 

Jc 

A = Y, S i^i < A ' 7) 

i* 1 

where 

{v\.}* ml i s an orthonormal set spanning the domain of A (A-8) 

{Ujjj.iis an orthonormal set spanning the range of A (A-9) 

{ s^} is a set of positive scalars (the singular 

values) (A-10) 

This is the SVD (Singular Value Decomposition) representation. 
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The special dyadic 


12 


* - £ 9k STi 

i =1 


(A- 11) 


constructed from an orthonormal set {c%} spanning the rnlimensional space is the identity dyadic in 
that space. Indeed, representing a in that frame as 

n 

<3=X) a i<7i (A- 12) 

i-l 


we get 


I-a 


u-i 


E \ E a j 9j 


V- 1 



(A- 13) 


Similarly, 


a -I = a 


(A-14) 


The properties established so far allow us to prove statement (2.5), namely, that the gradient 
of a position vector is the identity dyadic. To do that, we let the vector a be a position vector in the 
n-dimensional space. In this case, the gradient operator V is defined in terms of the components 
of a as 




i-i 


9i da d 


(A- 15) 
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Hence, 


Va = 




(A-16) 


But 



(A- 17) 


Hence, 


- n n 
Va = £ = £ gJi 

1, j = 1 - 1 * 1 


= I 


(A-18) 


and (2.5) is proved. 

The scalar product of two dyads is another dyad. Let 


F = cd 


(A-19) 


then 


D‘F = (ab) • (cd) 

= a (b -c) d 

= (jb-c) (ad) (A-20) 


Similarly, the product of two dyadics is another dyadic. 
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Finally, the inverse of a nonsingular dyadic N, denoted N 1 , satisfies 


AT 1 'N - I 


(A-21) 


N-N ~ 1 = I 


(A-22) 


We conclude with a derivation of the relationship between the dyadic representation of the 
operation of a linear operator and the more familiar matrix representation. Let the dyadic A operate 
on the vector q in then-dimensional input space (spanned by the orthonormal set { to yield 
the vector p in them-dimensional output space (spanned by the orthonormal set {e^ ” =1 ). The 
dyadic representation of this process is 


A-q - p (A-23) 

Now we replace q by I ~q , where X is the identity operator of the input space, and dot both sides of 
(A-23) by e s : 


( n 


e.-A 


E i i i i 


U-l 


•q = e ± -p 


(A-24) 


Extending the j-summation to cover all of the left side, we get 


E = (e^p) 

j=i 


(A-25) 


But this has essentially the matrix form we are after. Specifically, denoting 


(e^p) = p i = element 1 of the column matrix 
representing p 


(A-26) 
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(fj'Q) = Qj 


element j of the column matrix 
representing q 


(A-27) 


(e^A’fj) = A ij = element (i, j) of the rectangular matrix 
representing A 


(A- 28) 


we get the standard matrix representation of process (A-23) in the adopted orthonormal frames, 
namely, 

A ij^j = Pi (A ’ 29) 

(A-28) also provides the foundation for converting in the other direction. Let us use it to substitute 
forAjj in the following double summation: 


m n _ 


e i e i 


VJ-i 


U-i 


^ 5 


(A- 30) 


We have here two identity operators: the right one operates on the input space while the left one 
operates on the output space. In view of (A-7), (A-13), and (A-14), both of them may be dropped 
and we get 

A = JT (A-31) 

i^i J-1 

Note that, whereas the (A-31) representation of A is a sum of mn dyads (corresponding to 
the mn elements of the matrix A), its S VD representation is a sum of only k dyads (where k is the rank 
of A). This, however, requires knowledge of its singular value decomposition. 

A good treatment of the subject of dyadics is given in [9, chapter IV]. A more extensive 
treatment combined with the SVD formulation is given in [4], 
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APPENDIX B: THE FOCAL-PLANE DISTRIBUTION DUE TO AN OFF-BORESIGHT 

PLANE WAVE 


In Section III we dealt with the determination of the aperture-plane distribution from a given 
(measured) focal-plane distribution. Here we consider the reverse problem, namely, finding the 
focal-plane distribution due to an aperture-plane distribution generated by a (slightly) off-boresight 
incident plane wave. This is motivated by the need to find the rate of attenuation of the focal-plane 
field as we move away from the optical axis. This rate affects the selection of N— the number of 
rings in Fig. 2. 

A 

We start with the aperture-plane distribution u caused by a plane wave with propagation 
vector £ (see (2.1), (3. 9)) 

u = e i( *‘ f - ut > f B -l 


A 

Expressing this in terms of the normalized aperture-plane position vector R (see (3. 6), (3. 8)), we get 


U ( R ) = t) 


(B-2) 


Now we apply (3.13) to get the corresponding focal-plane distribution u (R) , namely. 


A A 

j j — g i + 1tR 2 -u> t) ^ ^ g 1 (L£ *R “2 it R *R) ^ ^ 


(B-3) 


where S refers to the aperture area and dS is an element of this area which we now express explicitly 
in its polar form 


dS = (rdi|f) dr = L 2 R dR d\ |r 


(B-4) 


Hence, 


A 



J?*0 ^ *0 


(B-5) 


A 

where R a is the normalized radius of the antenna aperture. 
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Next we turn our attention to the exponent in the integrand of (B-5). Substituting for Jc using 

A _ 

(2. 3), (2. 16), and noting that R \h 3 = 0 , we find that the exponent takes the following form: 

A 

-i2n [(L/X)n g + R]‘R (B-6) 


This suggests the introduction of the "offset position vector" b given by 

b = Fh g + f (B-7) 


where Fis the focal length of the antenna and r is the (unnormalized) position vector in the focal 
plane (3.10). The exponent (B-6) is now expressible in terms oiB, the normalized version 
of b , given by 


B = b/L = ( L/\)n g + R 


(B- 8) 


This simplifies (B-5), yielding 


U = 


2 n 


Jf 


L 2 

je=o t=o 


e - i2nB ' R R dR di|r 


(B-9) 


In evaluating this integral, we measure the polar angle ijr in reference to the vector B as shown 
in Fig. 18. This means that 


_ _ A 

B-R = BR cos ijr 


*5 

(B- 10) 


and (B-9) now takes the following form: 


U = 


A 

2 n 

L2 f 

f e -i2nBRcosy ^ 

A*' 

R= 0 

J 

_ 


A A 

R dR 


(B- 11) 
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The bracketed integral is now recognized as a multiple of the Bessel function of order zero ( J a ) (see 
formula 9. 1 .21 of [10]). Specifically, 


AAA 


u = 2uL y RdR 

2nBR M 

= [(L/B) 2 / (2w)] f xJ 0 (x) dx 


(B-12) 


But the integral appearing here is immediate and given by xJ x (x) where J x (x) is the Bessel 


function of order 1 (see formula 9.1.30 of [10]). This yields 


U = U(B) 


( 2 * f /) 


J 1 (2nBR a ) 

(2nBi? a ) 


(B-13) 


where r = is the unnormalized radius of the antenna aperture. 

a 

Note that for small x, 

cT 2 (x) * ( l/2 ) x ( ] x | < < 1) (B- 14) 


(see formula 9.1.7 of [10]). This implies 

lim J 1 (x)/x = 1/2 (B-15) 

X-0 

Furthermore, when we superpose the line y = (1/2) x on a plot of y = J x (x) (such as Fig. 

9.1 of [10]) we find that 

| ^(x) | £ | (l/2)x| (B- 16) 

We conclude that the function J t (x) /x as well as its magnitude assume their maximal values of 
1/2 at the origin. Applying this to (B-13), we get 

\CJ(B) | £ 17(0) = 7i r a 2 = area of aperture (B-17) 


Hence the normalized phrasing of (B-13): 


U(B) _ , 
17(0) 


J x (2nBR a ) 
(2nBR a ) 


(B- 18) 
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Note that ^differs from the field entity u only in phase (see (B-5)). Hence the maximal value 
of | u| occurs at B = 0 . Thus, the coordinate of the peak response r PEAK \s given by 


PEAK 


-Fn, 


(B-19) 


(see (B-7)) while the corresponding normalized coordinate (R PEAK ) is given by 

R peak ~ ~(B/X)n g 


(B-20) 


(see (B-8)). It follows that the value of u at this point ( u PEAK ) is given by 


U PEAK 


_ f/(0) [ 


(B-21) 


and finally 


u 

U PEAK 


= 2 


A 

J 1 {2nBR a ) e in{R'-UL/Un 3 ] 2 } 
(2nBR a ) 


(B-22) 


Note that the vectors R and B affect this expression only via their magnitudes. This means that the 
structure of the focal-plane field is describable in terms of two families of concentric circles (see 
Fig. 19); the constant-magnitude loci, which are circles of radius B centered at i? p£AJf and the 
constant-phase loci, which are circles of radius R centered at the origin. 

Applying formula 9.2.1 of [10] to (B-22), we find that for large values of the 
argument (2 tc BR a ) , 

I u/u PEAK \ < J871i/ (2nBR a ) 1 (2ttBR a > >1) (B-23) 

For the parameters of Fig. 7, this yields for the fourth ring 

\u/u PEAK \ < 0.007 {N = 4) (B-24) 

and thereby we are justified in ignoring terms of the fourth and higher rings. Note, however, that 
even when we use only one ring, the neglected terms are still quite small. Thus, 

l u / u PEA*l <0.021 (N = 2) (B- 25) 
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Before leaving this subject, it should be pointed out that (B-22) and Fig. 19 provide (at least 
in principle) an alternative approach to the DOA extraction. For example, if we subtract from the 
phase of each ( u/ u PEAK ) the known phase nR 2 and average these entities over all horns, we get an 
estimate of rig and hence the 0 coordinate of the source. Alternatively, one could ignore the phase 

A 

information and extract the parameter (2 nBR a ) from each focal-plane voltage. Each Bis the 
(normalized) distance from the known location of the measurement to the tip of the vector R PEAK 
(point C in Fig. 19). Obviously, if we succeed in locating point C, we have obtained R PEAK and, 
hence, ii s (B-20), and the problem is solved. It is interesting to note that, as in the adopted 
formulation of Section IV, we have here a purely geometric phrasing of the DOA extraction, namely, 
finding a point (C) from its distances (theB’s) to a given set of points. Here, we are dealing with 
distances to a point. There, we have been dealing with distances to a plane. 

An optimal combination of this approach with the one yielding rig is worth looking into. 
Whether the resulting algorithm would be practically feasible and advantageous is not currently 
known. 
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